{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# SciPy\n",
    "\n",
    "Examples taken from [Scipy Lecture Notes](http://www.scipy-lectures.org)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Loading default packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.display import display, Markdown\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.style.use(\"default\")\n",
    "#% matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# plot settings\n",
    "plt.rcParams['figure.figsize'] = [9.5, 4.25]\n",
    "plt.rcParams['figure.dpi'] = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from scipy import interpolate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Create sample signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "measured_time = np.linspace(0, 1, 10)\n",
    "noise = (np.random.random(10)*2 - 1) * 1e-1\n",
    "measures = np.sin(2 * np.pi * measured_time) + noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Perform `linear` and `cubic` interpolation with [`interp1d`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "linear_interp = interpolate.interp1d(measured_time, measures)\n",
    "interpolation_time = np.linspace(0, 1, 50)\n",
    "linear_results = linear_interp(interpolation_time)\n",
    "cubic_interp = interpolate.interp1d(measured_time, measures, kind='cubic')\n",
    "cubic_results = cubic_interp(interpolation_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Display results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(measured_time, measures, '.-', label='Measured signal')\n",
    "plt.plot(interpolation_time, linear_results, label='Linear interpolation')\n",
    "plt.plot(interpolation_time, cubic_results, label='Cubic interpolation')\n",
    "plt.xlabel('Time $t$ / s')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fast Fourier transforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from scipy import fftpack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Create a noisy sine signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = np.arange(0, 10, 0.01)\n",
    "u_t = (np.sin(2 * np.pi * 1 * t) + 0.5 * np.random.randn(t.size))\n",
    "plt.figure()\n",
    "plt.plot(t, u_t)\n",
    "plt.xlabel('Time $t$ / s')\n",
    "plt.ylabel('Voltage $u[t]$ / V')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Compute the Fast Fourier Transformation (real part - positive frequencies)\n",
    "\n",
    "The FFT $U[k]$ of length $N$ of the length-$N$ sequence $u[t]$ is\n",
    "defined as\n",
    "\n",
    "$$X[k] = \\sum_{n=0}^{N-1} e^{-2 \\pi j \\frac{k t}{N} } x[t]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "U_f = fftpack.rfft(u_t)\n",
    "U_f[0] /= u_t.size\n",
    "U_f[1:] /= u_t.size / 2\n",
    "f = fftpack.rfftfreq(u_t.size, 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(f, np.abs(U_f))\n",
    "plt.xlabel('Frequency $f$ / Hz')\n",
    "plt.ylabel('Absolute specturm |U(f)| / V')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimization and fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Create noisy data for a function $f(t) = a \\cdot \\sin\\left(b \\cdot t\\right)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x_data = np.linspace(-5, 5, num=50)\n",
    "y_data = 2.9 * np.sin(1.5 * x_data) + np.random.normal(size=50)\n",
    "plt.figure()\n",
    "plt.plot(x_data, y_data, '.', label='Measured data')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Create a prototype function for fitting and use [`curve_fit`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def test_func(x, a, b):\n",
    "    return a * np.sin(b * x)\n",
    "params, params_covariance = optimize.curve_fit(test_func, x_data, y_data, p0=[2, 2])\n",
    "display(Markdown(f'$a$ = {params[0]}, $b$ = {params[1]}'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(x_data, y_data, '.', label='Measured data')\n",
    "plt.plot(x_data, test_func(x_data, params[0], params[1]), label=r'Fitted curve $y(t) = a \\cdot \\sin(b \\cdot t)$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RL25",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
